/**
* @file MyMath.h
* @brief This file aims to deal with solving the matrix in Poisson problem.
*
* Firstly, a basic matrix solver is defined in this file, which is used to solve linear matrix problem with
* given sparse matrix, right-hand side of the equation. Secondly, derived solver class is designed to deal with
* building the required matrix and rhs with input physical parameters. Sparse matrix and related solving method are
* used in the code, for fast and accurate solution of Poisson equation.
*
* @author lunzhy
* @version 
* @date 2013-5-28   9:59
* @note This class is designed to provide all the math operation used in the process of solving Poisson problem
* @todo
*/
#pragma once
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/OrderingMethods>
#include <vector>
#include <iostream>
#include <vector>
#include "Utilities.h"
using namespace Utilities;


/// @brief Matrix solver provides basic method used in solving linear matrix equation
///
/// Sparse matrix is required in this class. SparseLU method for decomposition is used in solving linear matrix problem
class MatrixSolver
{
protected:
	Eigen::SparseMatrix<double> sparseMatrix;///< matrix A in Ax=b equation. It is filled in the derived class
	/// @brief SolveMatrix uses sparse matrix solver to solve linear equation
	/// 
	/// This method uses the SparseLU solver in Eigen to solve sparse matrix. The Map method in Eigen is used in this 
	/// method because the pass-by-reference parameters are std::<vector>. Eigen requirs its internal type of matrix
	/// vector. Map is used to deal with this problem.
	/// 
	/// @param std::vector<double> & rhs right-hand side of the equation
	/// @param std::vector<double> & solution the solution of the equation
	/// @pre
	/// @return int 0 for success, 1 for failure
	/// @note
	/// @see
	int SolveMatrix(std::vector<double> &rhs, std::vector<double> &solution);
};


/// @brief OneDimSolver is a derived class for solving one dimensional Poisson problem
///
/// The sparse matrix is build in this class. It can call the method in the base class to solve sparse matrix
class OneDimSolver : public MatrixSolver
{
public:
	OneDimSolver(int vertexNum);
	/// @brief SolveOneDim is used to solve specific one dimensional Poisson problem.
	/// 
	/// The process of solving one dimensional Poisson problem is organized in this method
	/// 
	/// @param std::vector<double> & rhs this parameter cannot be const probably because of the internal mechanism of Map is Eigen
	/// @param std::vector<double> & x
	/// @param const std::vector<double> & eps
	/// @param const std::vector<double> & deltaX
	/// @param const OneDimBC & left
	/// @param const OneDimBC & right
	/// @pre
	/// @return int
	/// @note
	/// @see
	int SolveOneDim(std::vector<double> &rhs, std::vector<double> &x, 
					const std::vector<double> &eps, const std::vector<double> &deltaX, const OneDimBC &left, const OneDimBC &right);
//////////////////////////////////////////////////////////////////////////
private:
	int oneDimSize;///< the number of the vertex in solving the problem
	/// @brief buildSparseMatrix is used to build specific sparse matrix with given physical parameters
	/// 
	/// This method fills the sparse matrix in the base class, which is used to solve the linear matrix equation.
	/// The sparse is refreshed in this method with given boundary conditions.
	/// 
	/// @param const std::vector<double> & eps the relative dielectric constant info of the mesh grid
	/// @param const std::vector<double> & deltaX mesh information
	/// @param const OneDimBC & left left boundary condition
	/// @param const OneDimBC & right right boundary condition
	/// @pre
	/// @return int
	/// @note the refreshing of rhs vector is done outside the solver class
	/// @see
	int buildSparseMatrix(const std::vector<double> &eps, const std::vector<double> &deltaX, const OneDimBC &left, const OneDimBC &right);
};


/// @brief TwoDimSolver is a derived class for solving two dimensional Poisson problem
///
/// The sparse matrix is build in this class. It can call the method in the base class to solve sparse matrix
class TwoDimSolver : public MatrixSolver
{
public:
	TwoDimSolver(int vertexNumX, int vertexNumY);
	/// @brief SolveTwoDim is used to solve specific two dimensional Poisson problem.
	/// 
	///
	/// 
	/// @param std::vector<double> & rhs cannot set &rhs to const, because it is used in SolveMatrix of base class, MatrixSolver
	/// @param std::vector<double> & x
	/// @param const std::vector<double> & meshX
	/// @param const std::vector<double> & meshY
	/// @param const My2DVector<double> & eps
	/// @param const TwoDimBC & left
	/// @param const TwoDimBC & right
	/// @param const TwoDimBC & up
	/// @param const TwoDimBC & down
	/// @pre
	/// @return int
	/// @note
	/// @see
	int SolveTwoDim(std::vector<double> &rhs, std::vector<double> &x,
					const std::vector<double> &meshX, const std::vector<double> &meshY, const My2DVector<double> &eps,
					const TwoDimBC &left, const TwoDimBC &right, const TwoDimBC &up, const TwoDimBC &down);
//////////////////////////////////////////////////////////////////////////
private:
	int twoDimSize;///< the total number of the vertex used in solving two dimensional problem
	int twoDimSizeX;///< the vertex number in X direction in terms of the two two dimensional distribution of the vertex
	int twoDimSizeY;///< the vertex number in Y direction in terms of the two two dimensional distribution of the vertex
	/// @brief buildSparseMatrix is used to build specific sparse matrix with given physical parameters
	/// 
	/// This method fills the sparse matrix in the base class, which is used to solve the linear matrix equation.
	/// The sparse matrix is filled regardless of the boundary conditions.
	/// 
	/// @param const std::vector<double> & meshX
	/// @param const std::vector<double> & meshY
	/// @param const My2DVector<double> & eps
	/// @pre
	/// @return int
	/// @note
	/// @see
	int buildSparseMatrix(const std::vector<double> &meshX, const std::vector<double> &meshY, const My2DVector<double> &eps);
	/// @brief refreshSparseMatrix refresh the sparse matrix with given boundary conditions
	/// 
	/// Currently, the boundary conditions of for entire boundary are considered.
	/// 
	/// @param const TwoDimBC & left
	/// @param const TwoDimBC & right
	/// @param const TwoDimBC & up
	/// @param const TwoDimBC & down
	/// @pre
	/// @return int
	/// @note Partial boundary condition is important when the contact is not set on the entire boundary.
	/// @see
	int refreshSparseMatrix(const TwoDimBC &left, const TwoDimBC &right, const TwoDimBC &up, const TwoDimBC &down);
};